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CHAPTER  1 
INTRODUCTION 

Automatic  error  control  has  been  studied  and  implemented 
since  the  late  fifties  (15) •      Over  the  past  decade  refinements 
have  been  made  and  further  analytic  results  have  come  to  light 
(^s  5 »  8,  9,  16).   It  appears  that  the  state  of  the  art  suggests 
that  designers  and  implementers  should  try  to  embed  these  rules 
into  their  software  systems.   This  provided  the  motivation  for 
this  thesis,  which  was  to  implement  the  rules  of  "unnormalized 
significant  digit  arithmetic"  of  Metropolis  (15>  9)  and  Ashen- 
hurst  (2,  4)  and  embed  them  in  the  array  language  OL/2  (20). 
However,  these  software  packages  are  also  adaptable  to  other  com- 
patible software  systems,  since  they  are  written  in  assembly 
language  for  the  IBM  360/370  machines. 

Numerical  errors  in  digital  computation  arise  from  three 
principle  sources:   "inherent"  errors  which  result  from  inaccur- 
ate data,  "analytic"  errors  which  result  from  replacing  an  in- 
finite process  with  a  finite  process,  and  "generated"  errors 
which  result  from  using  finite-precision  arithmetic  (1).   Often, 
computed  results  are  reported  without  specifying  their  accuracy 
because  inherent  and  generated  errors  have  been  treated  rather 
lightly.   Presently,  most  floating-point  arithmetic  is  performed 
in  normalized  form  which  leads  to  the  retention  of  nonsignificant 
digits  in  final  computed  results.   There  is  no  control  of 


significant  digits  when  experimental  data  is  used  as  input. 
Even  with  precise  input,  ill-conditioning  or  instability  in 
the  programmed  algorithm  may  cause  erosion  of  precision,  thus 
a  specification  of  error  or  significance  should  accompany  out- 
put numbers. 

As  a  consequence  of  its  practical  importance,  analytical 
investigations  into,  as  well  as  implementation  of,  significance 
control  has  been  made  by  several  earlier  investigators  on  sev- 
eral different  machines.   Notable  among  these  are  Gardiner  and 
Metropolis  (10),  Goldstein  (11),  Ashenhurst  (  3  ),  and  Bright' 
(  7  ) .   Each  agree  to  the  two  rules  of  significant  digit  arith- 
metic:  (1)  for  addition  and  subtraction,  the  resulting  exponent 
is  equal  to  the  larger  of  the  input  exponents,  and  (2)  for  mul- 
tiplication and  division,  the  number  of  significant  digits  in 
the  result  is  equal  to  that  of  the  less  significant  operand. 
They  differ,  however,  in  their  implementations  of  the  internal 
representation  of  significant  digit  numbers.   Because  each  mach- 
ine employed  a  different  exponent  base,  corresponding  significant 
adjustment  rules  for  library  functions  also  reflect  this  differ- 
ence . 

All  numerical  quantities  may  be  regarded  as  either  exact 
or  approximate.   However,  the  finite  presentation  of  number  in 
computer  requires  a  further  characterization.   If  an  exact  number 
can  be  represented  in  some  preassigned  number  of  digits  (the 
computer  word  length),  then  it  is  said  to  be  precise .   Other 
exact  numbers  belong  to  the  approximate  class  are  called 


3 
imprecise .   Using  these  two  categories  and  base  of  two,  Gardi- 
ner and  Metropolis  (  9  )  introduced  the  ABC  (analyzed  binary 
computation)  form  of  arithmetic.   The  goals  of  ABC  arithmetic 
are:   (1)  to  determine  if  the  computed  result  is  precise  or  un- 
precise,  and  (2)  to  specify  the  number  of  significant  digits  in 
the  imprecise  result.   Thus,  ABC  arithmetic  uses  the  rules  of 
significant  digit  arithmetic  and  defines  an  additional  set  of 
rules  for  combining  both  precise  and  imprecise  numbers.   When- 
ever we  use  the  term  "significant  digit  arithmetic"  we  mean  to 
include  the  new  rules  of  ABC  arithmetic  as  well.   The  present 
work  follows  the  approach  of  Gardiner  and  Metropolis  (  9  )  >  ex- 
cept that  the  base  two  rules  must  be  replaced  by  base  sixteen 
rules,  since  the  implementation  is  for  the  IBM  360/370  machines. 

Overview  of  the  Thesis 

The  next  chapter  defines  the  basic  rules  for  combining 
precise  and  imprecise  numbers  in  the  context  of  ABC  arithmetic. 
Two  kinds  of  zeroes,  true  and  relative,  are  also  discussed  in 
this  context,  along  with  the  internal  format  for  imprecise  and 
precise  numbers. 

Chapter  3  discusses  the  significant  digit  adjustment  rules 
for  elementary  functions.   In  particular,  Richtmyer's  adjustment 
rule  for  functions  is  derived.   The  discussion  follows  the  rules 
that  were  prescribed  by  Ashenhurst  (4),  with  suitable  changes 
for  base  sixteen  rather  than  base  two  numbers. 

The  fourth  chapter  explains  the  format  of  input/output  for 
precise  and  imprecise  numbers.   The  significance  number  base 
conversion  algorithms  of  Kanner  (13)  will  be  presented. 


Chapter  5  briefly  explores  some  extensions  to  the  present 
work,  and  also  some  of  the  problems  which  the  uninitiated  user 
must  be  aware  of  when  using  significance  arithmetic . 


CHAPTER  2 
PRECISE/IMPRECISE  ARITHMETIC 

2.1  Unnormalized  Floating  Point  Numbers 

A  floating  point  number  may  be  represented  by  a  pair 
(e,  f),  where  e  is  the  integral  exponent  and  f  the  real  valued 
fraction  part.   Most  of  today's  computers  use  normalized  arith- 
metic in  which  the  floating  point  number  has  its  non-zero  frac- 
tional part  satisfying  r   ^_  f  <  r  ,  where  r  is  the  number  base. 
This  format  alone  cannot  specify  the  number  of  significant  digits 
in  the  fraction.  By  unnormalizing  the  number,  it  is  possible  to 
convey  the  significant  digits  in  the  fraction,  since  (e  +  m, 
r   f)  represents  the  same  number  as  (e,  f)  for  any  integer  m. 
The  integer  m  describes  the  number  of  leading  zeroes  in  the 
fraction.   By  definition,  therefore,  imprecise  numbers  are  re- 
presented internally  as  unnormalized  numbers  with  at  least  one 
leading  zero.   Precise  numbers,  on  the  other  hand,  are  repre- 
sented internally  as  normalized  numbers.   This  distinction  be- 
tween precise  and  imprecise  numbers  forms  the  basis  for  the  rules 
of  ABC  arithmetic  (9)  which  are  discussed  subsequently. 

2.2  True  Zero  and  Relative  Zero 

In  significant  digit  arithmetic  there  are  two  kinds  of 
zeroes:   relative  zero,  which  has  no  significant  digits  and  true 
zero,  which  has  an  infinite  number.   Both  zeroes  have  zero  frac- 
tional parts,  and  obviously,  there  is  no  sign  attached  to  the 


6 
fractions.   True  zero  is  designated  by  using  the  smallest  ma- 
chine representable  exponent,  in  general,  this  is  just  the  usual 
machine  representation  of  zero.   Relative  zero,  on  the  other 
hand,  may  have  any  exponent  value.   Clearly,  the  rules  for  operat- 
ing with  these  two  zeroes  is  quite  different.   True  zero  is  a 
precise  number  in  the  sense  that  it  acts  like  zero.   In  contrast, 
relative  zero  is  a  degenerate  imprecise  number  with  no  signifi- 
cant digits  but  with  an  exponent  field  that  indicates  an  "order 
of  magnitude."   Since  relative  zero  contains  some  information 
about  a  number  besides  complete  loss  of  significance,  it  seems 
appropriate  to  consider  separate  rules  for  operating  with  it . 
It  should  be  clear  that  a  relative  zero  may  occur  during  a 
computational  process  as  the  result  of  very  basic  operations. 
For  instance,  a  relative  zero  may  occur  when  two  nearly  equal 
imprecise  numbers  are  subtracted.   If  the  original  operands 
had  a  large  exponent,  then  the  relative  zero  has  a  large  ex- 
ponent.  This  indicates  that  a  relative  zero  is  not  necessarily 
close  to  zero,  but  is  merely  a  bound  on  the  result. 
2.3   Internal  Representation  of  Numbers 

The  present  work  was  implemented  on  a  IBM  360/75  computer 
using  long  floating  point  numbers.   These  numbers  use  64  bits: 
7  bits  for  the  base  16  exponent  and  57  bits  for  the  fraction. 
(See  Figure  1.)   The  zero-th  bit  is  used  for  the  sign  of  the 
fraction.   The  exponent  is  represented  in  excess-64  code  with 
limits  of  -63  and  +  64.   True  zero  has  an  exponent  of  -64  and 
has  zeroes  in  all  64  bits.   (This  Is  the  standard  floating  point 


7 
zero.)   Relative  zero  is  represented  by  a  zero  fraction,  a  zero 
in  the  sign  position,  and  arbitrary  exponent  e  in  the  range 
-63  <_   e  £  63.   The  sign  of  relative  zero  is  neither  positive  or 
negative . 

A  precise  number  is  represented  as  normalized  floating 
point  number  with  14  hexadecimal  digits  in  the  fractional  part. 
Only  precise  numbers  are  normalized.   The  range  of  IBM's  long 

7  7         -t7  C\ 

floating  point  number  is  from  10    to  10    . 

Since  an  imprecise  number  is  an  approximate  quantity,  it 
is  represented  with  fewer  fractional  digits.   The  first  4-bit 
fractional  byte  is  always  zero  and  the  last  hexadecimal  digit 
is  always  defined  to  be  a  nonsignificant  guard  byte.  Thus,  the 
imprecise  fractional  field  uses  the  middle  12  hexadecimal  digits. 

Because  imprecise  numbers  are  unnormalized,  they  can  only  repre- 

-77       +49 
sent  numbers  between  10    and  10    .   The  guard  byte  is  used 

in  input/output  conversion  and  also  serves  as  a  guard  against 
round-off  errors.   Kanner  (13)  points  out  that  if  input  conver- 
sion followed  by  output  reconversion  is  to  be  an  identity  opera- 
tion at  least  one  guard  bit  must  be  kept.   Because  the  exponent 
base  is  sixteen,  it  is  necessary  in  this  implementation  to  use 
an  entire  byte  of  4  bits. 
2.4   Base  Operations 

In  this  section  the  operations  of  addition  (including  sub- 
traction), multiplication,  and  division  in  ABC  arithmetic  are 
described.   If  one  of  the  operands  is  imprecise,  then  the  result 
is  imprecise  or  a  relative  zero.   The  significant  digit  rules  are 
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fraction  field  (8:63) 
exponent  field  (1:7) 
fraction  sign  bit  (0) 


a)   IBM  360/370  long  floating  point  format 


1 

1 

guard  byte  (60:63) 


.leading  zero  byte  (8:11) 


b)   Imprecise  number 


Figure  1 .   Internal  number  format 


9 
used  for  imprecise  numbers.   Let  X,  =  (e,,  f -,  )  and  Xp  =  (e^,  fp) 
be  the  operands,  and  X^  =  (e~,  f-O  the  results,  then  the  rules 
of  significant  digit   arithmetic  are  defined  as  follows: 

1.  For  addition  or  subtraction  e_  =  max  {e, ,  ep}, 

2.  For  multiplication  or  division  m.3  =  max  {iru  ,  nip} 
where  m  is  the  number  of  leading  zeroes  in  the  fraction.   If  a 
zero  fraction  results,  the  exponent  of  the  relative  zero  is 
fixed  according  to  the  rules  of  regular  arithmetic.   That  is, 

in  addition  e~  is  the  maximum  exponent  of  operands,  in  multipli- 
cation e,  is  the  algebraic  sum  of  the  exponents,  and  in  division 
e~  is  the  difference  between  the  exponent  of  the  numerator  and 
the  exponent  of  the  denominator.   There  are  some  obvious  prob- 
lems with  using  the  360/370  floating  point  arithmetic.   If  the 
unnormalization  of  imprecise  numbers  Is  not  possible,  then  the 

ABC  arithmetic  routines  signal  exponent  overflow.   This  is  why 

+4  9 
imprecise  numbers  have  an  upper  limit  of  10    instead  of  the 

+75 
10    of  the  precise  number. 

A  necessary  but  not  sufficient  condition  for  a  precise  re- 
sult is  that  both  operands  are  precise.   The  result  would  be 
imprecise  if  it  cannot  be  represented  exactly  in  the  long  float- 
ing point  word.   For  multiplication  and  division,  the  low  order 
part  of  the  immediate  product  is  generated  and  examined  for  any 
non-zero  bits.   (The  360/370  machines  generate  only  the  higher 
order  part.)   Since  the  product  is  almost  twice  the  length  of 
input  operand,  a  quick  test  for  a  precise  multiply  is  to  check 
if  both  precise  inputs  are  contained  in  the  short  floating  point 


10 
number  format  which  has  32  bits  and  6  hexadecimal  fractional 
digits.   For  addition,  the  precise  check  is  to  subtract  the 
result  by  the  operand  with  the  larger  exponent.   This  subtrac- 
tion would  have  no  shifting  of  fraction  to  align  the  exponents 
The  result  is  then  compared  with  the  other  input  operand.   If 
they  are  not  equal,  fractional  digits  were  lost  in  the  addition; 
so  the  result  is  imprecise.   The  special  case  of  a  true  zero 
resulting  from  the  subtraction  of  the  precise  number  by  itself 
is  checked  by  comparing  the  input  fractional  digits  for  equality. 

As  noted  earlier,  true  zero  is  the  zero  element  in  ABC 
arithmetic.   In  addition,  it  is  the  identity  element.   A  true 
zero  product  results  if  one  of  the  operands  is  a  true  zero. 
Division  by  true  zero  is  signaled  as  a  divide  exception  and  a 
true  zero  dividend  results  in  a  true  zero  quotient. 

On  the  other  hand,  the  rules  for  operating  with  relative 
zero  are  very  different.   If  both  operands  are  relative  zeroes, 
then  the  exponent  of  the  resulting  relative  zero  exponent  must 
be  adjusted.   In  extreme  cases,  an  exponent  overflow  in  a  rela- 
tive zero  raises  an  error  condition  while  an  exponent  underflow 
is  replaced  with  the  smallest  relative  zero.   In  multiplication 
and  division,  the  result  is  generally  a  relative  zero  if  one  of 
the  operands  is  a  relative  zero.   It  is  interesting  to  note  that 
the  quotient  of  a  non-zero  divided  by  a  relative  zero  is  defined 
as  a  relative  zero.   In  this  case,  multiplication  of  quotient 
and  divisor  does  not  give  the  dividend. 
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Table  1 
ABC  Multiplication 


°T 

°R 

XI 

xp 

°T 

°T 

°T 

°T 

°T 

°R 

°T 

°R 

°r 

°R 

XI 

°T 

°R 

°R:XI 

°R:XI 

xp 

°T 

°R 

°R:XI 

Xj:Xp 

0™  =  true  zero       XT  =  imprecise  non-zero 
0R  =  relative  zero   Xp  =  precise  non-zero 


Table  2 
ABC  Division 

°T 

°R 

xi 

Xp 

°T 

00 

00 

oo 

oo 

°R 

oT 

°r 

°r 

°r 

XI 

oT 

°r 

°r:Xi 

°R:XI 

XP 

oT 

°r 

°R:XI 

Xx:Xp 

Orp  =  true  zero       XT  =  imprecise  non-zero 
0R  =  relative  zero   Xp  =  precise  non-zero 
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The  addition  of  a  relative  zero  and  non-zero  may  have 
either  an  imprecise  result  or  a  relative  zero.   If  the  exponent 
of  the  non-zero  operand  is  greater  than  or  equal  to  the  exponent 
of  the  relative  zero,  then  the  sum  is  imprecise.   In  particular, 
if  the  non-zero  operand  is  precise,  then  the  sum  Is  the  rounded 
imprecise  representation  of  the  precise  non-zero  operand; if  the 
non-zero  operand  is  Imprecise,  then  the  sum  is  the  imprecise 
non-zero  operand  adjusted  (signifying  loss  of  significant  digits) 
by  the  difference  of  exponents  of  the  non-zero  operand  and  rela- 
tive zero.   Clearly,  if  the  exponent  of  the  relative  zero  is 
greater  than  the  exponent  of  the  non-zero  operand,  then  the  sum 
is  just  the  relative  zero. 

Table  3 
ABC  Addition  and  Subtraction 


°T 

°R 

XI 

XP 

°T 

°T 

°R 

XI 

XP 

°R 

°R 

°R 

°R:XI 

°R:XI 

XI 

XI 

°R:XI 

°R:XI 

°R:XI 

xp 

xp 

°R:XI 

°R:XI 

Om  :  X-j- :  Xp 

0m  =  true  zero 


XT  =  imprecise  non-zer< 


0R  =  relative  zero     Xp  =  precise  non-zero 
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The  rules  just  described  for  significant  digit  are  arith- 
metic implemented  in  assembly  language  as  subroutines.   The 
entry  point  names  are  @  0L2ADD,  @  0L2SUB,  PIMPY,  and  PIDIV. 
(See  Appendix  B  for  the  program  listings.)   Each  routine  employs 
a  round  on  the  digit  to  shift  off  due  to  unnormalization  of  the 
result  to  the  correct  number  of  significant  digits. 

There  are  some  problems  with  significant  digit  arithmetic 
using  unnormalized  numbers.   One  aspect  will  be  discussed  here 
and  others  will  be  treated  in  Chapter  5.   The  significance  rule 
for  multiplication  and  division  calls   for  retention  of  the 
same  number  of  significant  digits  as  in  the  operand  with  the 
fewest  number  of  significant  digits.   This  rule  can  result  in 
the  loss  of  significant  digits.   For  instance,  consider  the 
following  multiplication:     01  x  09  =  09  and  09  *  09  =  81. 
In  both  cases,  the  rule  specifies  one  significant  digit.   For 
the  first  case,  it  is  correct;  however,  in  the  second  case  the 
second  digit  "1"  would  be  discarded.   Thus,  a  corrective  shift 
or  an  auxiliary  rule  is  necessary  for  multiplication  and  divi- 
sion.  The  auxiliary  rule  for  multiplication  is:   if  e~  =  e,  +  e?, 
then  itu  =  m~  -  1,  where  e.'s  are  the  exponents  of  the  normalized 
operands.   This  corrective  right  shift  for  mutliplication  re- 
sults in  gaining  one  significant  digit.   Since  division  is  the 
inverse  operation  of  multiplication,  the  division  auxiliary  rule 
is:   if  e.  =  e,  -  e-,  then  m~  =  m_  +  1 .   Here  the  corrective 
left  shift  in  division  results  in  losing  a  significant  digit. 


in 

Routines  @  0L2MPY  and  @  0L2DIV  are  Identical  to  PIMPY  and 

PIDIV  except  for  the  incorporation  of  the  above  auxiliary  rules. 

The  next  chapter  will  examine  the  rules  of  significant 
digit  arithmetic  for  elementary  functions  using  base  sixteen. 
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CHAPTER  3 
ADJUSTMENT  RULES  FOR  FUNCTIONS 

3.1   Function  Error  Propagation 

Let  F  be  a  function  of  one  argument  and  let  (e-^,  f ,  )  be  the 
argument  value  at  which  F  is  to  be  evaluated.   Clearly,  this 
requires  the  specification  of  two  functions,  g  and  h: 

f2  =  g(ei,  f1) 

e.2  =  h(e1,  f1) 

e_        e^ 
such  that  r  ^f?  =  F(r   ,f-,).   This  relation,  however,  does  not 

imply  anything  concerning  the  number  of  significant  digits  In 
the  result  (e?,  fp)  since  (e?  +  m,  r   fp)  also  satisfies  it. 
The  adjustment  criteria  must  be  based  on  further  considerations. 
Given  a  number  (e-,  ,  f -,  )  ,  let  f-.  be  regarded  as  an  approxi- 
mation to  some  true  value  f-i^.  The  error  in  the  function  of  the 

T 
argument  is  defined  by  6,  =  f,  -  f,  ,  while  the  error  in  frac- 
tion of  the  function  is  defined  by  &~   =  g(e,  ,  f -,  )  -  g(e-.  ,  f-.  )  . 
If  we  use  the  generalized  mean  value  theorem,  then  we  have  the 
estimate  |6p|  =  a|6, |.   The  error  propagation  associated  with 
a  given  g  is  thus  expressible  in  the  terms  of  the  "amplification 
factor"  a.   This  is  the  procedure  which  is  described  by  Ashen- 
hurst  ( 4 )  . 

Adjustment  of  the  function  value  (ep,  fp)  is  made  so  that 
a  is  within  a  factor  of  16  (since  we  are  using  base  sixteen)  of 
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of  any  predetermined  value.   If  a  Is  approximately  one,  then 
an  adjustment  is  said  to  be  in  accordance  with  a  "significant 
digit"  criterion.   Notice  that  this  does  not  imply  that  all  the 
computed  digits  of  f„  are  actually  significant.   Rather,  it 
means  that  the  error  in  f«  due  to  the  error  in  f-,  is  of  the 
same  order  of  magnitude  in  terms  of  the  adjustment  used.   Both 
f..  and  fp  contain  roughly  the  same   number  of  "non-significant" 
digits . 

3.2   Propagation  of  Relative  Error:   Richtmyer  Adjustment  Rule 

T 
Suppose  x  is  a  number  approximating  a  true  value  x   not 

equal  to  zero.   The  relative  error  between  the  approximate  and 

T    T 
true  values  may  be  defined  as  r(x)  =  (x  -  x  )/x  .   Similarly, 

the  relative  error  in  the  function  F  may  be  defined  as  r„(x) 

=  (F(x)  -  (F(xT))/F(xT) .   If  F(x)  is  dif f erentiable  at  x,  then 

T       T 
r„(x)  =  (F'(x)  (x  -  x  ))/F(x  ).   Using  the  definition  of  r(x), 

T  T  T 

we  have  rp(x)  =  (x  F ' (x)r (x) )/F(x  ).   If  (x  -  x  )  is  small,  we 

T 
replace  x   with  x  to  get: 

,  X     x  P'  (x)    ,  s 

r-o(x)  =  — -p /  \      r(x). 
■pv  /      F(x) 

By  definition  l6"m-1  <  |f|  <  l6"m,  thus  it  follows  that 
l6m|6|  <  \^\    <  l6m   |6|.   Since  this  inequality  may  be  applied 
to  both  the  argument  and  function  values,  using  respectively 
f .,  ,  m,  ,   .  ,  and  fp,  nip,   ~ ,  we  can  write 

l6m2-m1-l   6. 


62fl   <  i6m2~ml+1  |62 


6lf2 
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Thus  , 


62fl 


6lf2 


I6e2fi.    I6eif. 


I6e2f 


F(x)  -  F(x  ) 


I6el6. 


T 


FTry 


T 


x  -  x 


rp(x) 


provided  that  6,  and  6?  are  sufficiently  small.   An  approximate 
inequality  which  holds  to  a  first  order  approximation  may  be 
formed : 

x  F' (x) 


j|  <  l6ml-m2 


which  is  equivalent  to 


TJJJ 


<   16a 


l6m1-m2-l 


x  F' (x) 
x 


<  a  <  I6ml-m2 


-mo+1 


x  F'(x) 


The  above  inequality  bounds  a  as  a  function  of  the  argument  x 
and  the  adjustment  (mp-m-,).   Taking  the  base  16  logarithm  of 
all  members  of  the  inequality  (this  is  permissible  since  log, r 
is  a  monotonic  increasing  function)  gives 


log 


16 


x  F' (x) 
TXT) 


-  Am-l  <  loi 


a  <  log- 


x  F  T  ( x ) 
TZx) 


-  Am+1 . 


'16  "  "  ^616 

With  a  "significant  digit"  criterion  a=l  the  above  inequality 
can  be  written  as 


-1  <  lo 


gl6 


x  FT (x) 


-  Am  <  1 
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Thus,  an  approximate  adjustment  criterion  Is 

x  F» (x) 


Am  *  l°gl6 


TT7T 


R.  D.  Rlchtmyer  had  suggested  the  above  adjustment  rule  for 
unnormalized  arithmetic  as  well  as  for  functions.   The  above 
deviation  is  identical  to  the  one  of  Ashenhurst  (  l]   )  except  base 
2  was  replaced  by  base  16. 

Notice  that  in  actuality  Am  is  between  -1  and  +1  (base  16), 
a  defect  which  Ashenhurst  points  out  can  be  troublesome  since 
it  implies  that  16    <  a  <  16.   He  shows,  however,  that  if  it 
is  assumed  that  16  If,  is  uniformly  distributed  in  the  interval 
(16   ,  16)  the  expected  value  of  a  is  unity.   Therefore,  on  the 
average,  Richtmyer's  Rule  provides  the  optimal  adjustment. 
3.3   Application  of  Richtmyer's  Adjustment  Rule 

This  section  investigates  the  adjustment  rules  for  various 
functions  in  base  16 .   The  difference  between  corresponding 
function  adjustment  rules  for  base  16  and  base  2  will  be  explained. 

Raising  to  a  power,  F(x)  =  x  ,  where  n  is  any  number  gives: 


Am  =  lo 


sl6 


/   n-lv 
x(nx    ) 


n 


=  logl6  |n 


Since  for  the  most  practical  purposes  [n|  «  16,  an  approximate 
adjustment  rule  is  Am  =  0.   For  n=l  and  n=2,  the  adjustment  rule 
is  m  =m,  ,  which  means  that  the  number  of  leading  zeroes  in  the 
normalized  fraction  remains  unchanged.   For  the  square  root 
function  where  n=l/2,  the  adjustment  rule  must  be  also  mp=m,  to 
conform  to  that  of  square  function.   However,  in  base  2, 
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Am  =  logp  |n|  dictates  m2=m, -1  for  n=l/2  and  nip=in,+l  for  n=2. 
Therefore,  the  adjustment  in  base  16  is  relatively  coarse  com- 
pared to  base  2.   That  is,  in  base  2  we  would  actually  adjust 
m  for  small  values  of  n,  whereas  in  base  16  the  adjustment  would 
be  unnecessary. 

For  exponentials  and  logarithms,  F(x)  =  a  =e      , 
Richtmyer's  Rule  specifies 


Am  =  logl6 


x 


x(a  In  a) 


=  logl6|x|+  logl6  | In  a| 

Recalling  16"  ■*-"   <_  |  f  -,  |  <16~  ■*■ ,  the  adjustment  Am  may  be  ex- 
pressed as  either  Am  =  (e-.-m-.-l)  +  log-./-  |  In  a  |  or  as 
Am  =  (e-,-m-,)  +  log-,/-  |  In  a|.   Once  again,  if  a  uniform  distribu- 
tion of  16  if.,  is  assumed,  these  two  equations  may  be  averaged 
to  give  Am  =  e,-m..-l/2  +  log-,/-  |  In  a|.   Thus,  the  adjustment  rule 
for  an  arbitrary  a  is  m?  =  e,  -  1/2  +  log-,/-  |  In  a|  .   Clearly, 
the  right  hand  side  of  the  equation  must  be  rounded  to  an  in- 
teger.  In  a  similar  fashion,  the  logarithm  function  F(x) 

=  log0  (x)  where  x>0  should  be  adjusted  by  the  rule: 
a 

m2  =  el  +  "'"/^  "  l°Si  6  I  ^n  al*   Notice  that  the  only  difference 
of  adjustment  rules  for  exponentials  and  logarithms  is  the 
factor  of  one-half. 

More  complicated  functions  such  as  trigonometric  functions, 
hyperbolic  functions,  and  the  error  function  may  be  expressed 
as  a  power  series  in  x .   As  examples,  consider  the  following 
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expansions   with    0   <_  x   <_  tt/^I  : 

2  i\ 
x           x 

cos    x    =    l-py+jp-+.     .     . 

3  5 

X  X 

sin  x   =   x   -   =-r  +  -p-p  +    .     .     . 

j  !  b  I 

If  we  assume  as  Ashenhurst  does  that  the  error  in  the  function 
can  be  approximated  by  the  error  in  the  first  term  of  the  series, 

the  sine  can  be  thought  of  as  a  power  series  in  x  and  the  cosine 

2 
can  be  thought  of  as  a  power  series  in  x  .   Using  the  previous 

derived  rule  for  P(x)  =  x  ,  namely, 


rru  =  iru  +  log-,/-  |n| 

we  have  for  sine  (n=l)  nip=m,  and  for  cosine  (n=2) 
nip  =  m,  +  log-,/-  2.   The  last  adjustment  rule  can  be  approximated 
by  mp=m-,  .   Ashenhurst  derived  the  cosine  adjustment  in  base  2 
as  nip  =  2e-.  -  m,  .   Again,  this  difference  is  due  to  the  different 
bases  used.   Using  the  first  term  for  an  approximation,  the  ad- 
justment rule  mp=m-,  also  holds  for  hyperbolic  functions,  expon- 
tial  functions,  and  the  error  functions. 

Significant  digit  function  routines  have  been  implemented 
as  subroutines  which  call  the  corresponding  IBM  PL/1  library 
functions.   Ashenhurst  in  ( 4  )  suggested  this  strategy  of 
carrying  out  the  function  computation  by  normalized   and  then 
adjusting  the  result  at  the  end.   The  adjustment  rules  for 
imprecise  arguments  for  functions  in  base  16  falls  into  two 
equations — either  mp=m-,  or  m?=  |e-,  |  . 


Table  4 
ABC  Function  Rules 
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Function 


Base  2 
Adjustment 


Base  16 
Adjustment 


Sin 

Cos 

Tan 

Square  Root 

Exponent :  e 


x 


Log 
Log, 


10 


'2 
Sinh 

Cosh 

Tanh 

Arctan 

Error  Function 

Complement 

Error  Function 


mp=m, 

m„=  I  2e-,  -in-, 
nip=m, 
mp=m, -1 

m2= I el-1 1 


m2-|e1| 

m?=  |  e-,  +1 


mp  I  ei 


mp=m, 

m?=  |  2e-,  -m-,  | 


mp=m-| 


m2=ml 


m2=ml 


m„=in. 


mp=m, 
mp=m. 


mp=m.| 


m0=m. 


m2=|e1 


mp=  I  e-i  I 
m2=|e1| 


nip=m. 


mp=m-| 


mp=m-. 


mp=iru 


mp=m-. 


mp=m. 


The  previous  discussion  assumed  the  arguments  were  im- 
precise or  true  zero.   The  subroutines  (except  for  the  square 
root)  always  return  an  imprecise  number  even  when  the  argument 
is  precise.   The  justification  for  this  action  is  that  values 
of  the  function  can  only  rarely  be  represented  as  a  precise, 


22 

number,  for  Instance,  sin  30°  =  .5  and  cos  0°  =  1 .   For  a  rela- 
tive zero  argument,  the  argument  exponent  is  examined.   If  the 
exponent  is  positive,  the  same  relative  zero  is  returned  for 
the  function  value.   Relative  zeroes  with  non-positive  exponents 
are  treated  as  if  the  argument  is  near  zero.   In  particular, 
an  imprecise  one  with  one  significant  digit  is  returned  for  the 
exponential,  cos,  cosh,  and  error  functions,  while  a  relative 
zero  with  exponent  of  +14  is  returned  for  the  sine,  tan,  sinh, 
tanh,  arctan,  and  complement  error  functions.   An  error  condi- 
tion is  signaled  when  relative  zeroes  with  non-positive  exponents 
are  arguments  to  the  log  functions.   For  the  square  root  func- 
tion, the  same  relative  zero  is  returned  as  the  functional 
value.   This  implementation  of  relative  zero  as  an  argument  for 
the  above  functions  does  not  account  for  repeated  functional 
evaluations  with  relative  zero.   Thus,  the  treatment  of  relative 
zero  should  be  studied  further  and  refined. 

This  completes  the  discussion  of  implemented  significant 
digit  adjustment  rules  for  functions.   The  next  chapter  is  con- 
cerned with  the  input  and  output  convention  for  precise  and 
imprecise  numbers. 
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CHAPTER  4 
INPUT  AND  OUTPUT 

This  chapter  will  deal  with  the  problems  of  input  and 
output.   We  shall  also  discuss  algorithms  for  number  base  con- 
version as  described  by  Kanner  (13)  for  significant  digit 
arithmetic . 
4.1   Objectives  of  Significant  Digit  I/O 

The  main  objective  of  I/O  for  significant  digit  arithme- 
tic is  to  obtain  a  faithful  representation  of  significance 
under  the  operation  of  base  conversion  from  binary  to  decimal, 
as  well  as  the  inverse  of  this  operation.   For  instance,  the 
number  2.00  on  input  should  be  immediately  converted  on  output 
(without  any  intervening  computation)  as  2.00,  and  not,  say, 
as  1.99*   Kanner  points  out  this  can  be  achieved  only  if  at 
least  one  guard  bit  is  maintained  in  the  floating  point  frac- 
tion.  In  this  implementation,  as  already  pointed  out,  a  guard 
byte  (4  bits)  is  always  available.   This  multipurpose  guard 
byte,  whose  bits  are  considered  to  be  nonsignificant,  also  serves 
to  separate  accumulated  rounding  error  from  the  error  inherent 
to  the  computation.   Since  I/O  is  part  of  the  human-computer 
interfacing  problem,  there  should  be  ample  diagnostic  aids  to 
help  the  user  interpret  the  computer  outputs.   The  I/O  routines 
of  the  present  work  were  written  with  the  above  objectives  in 
mind. 
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4.2   Decimal  to  Binary  Conversion 

The  following  algorithms  of  Kanner  (13)  are  based  upon 
the  strategy  that  correct  significance  can  be  retained  in  the 
conversion  of  integers.   Briefly,  the  procedure  involves  multi- 
plying the  number  by  that  power  of  ten  which  is  necessary  to 
express  the  significant  portion  as  an  integer,  then  this  inte- 
ger is  converted  to  its  intermediate  internal  representation  in 
the  desired  number  system,  and  finally,  this  result  is  divided 
by  the  same  power  of  ten. 

The  following  discussions  have  been  adapted  to  base  16 

instead  of  base  2  as  discussed  in  Kanner  (13).   The  binary 

representation  of  N  can  be  obtained  from  the  decimal  integer 

representation  w,Wp...w  ,w  by  considering  the  following 

scheme:   a  =0,  a.=10a.  -,    +  w.  ,  N=a  .   The  long  floating  point 
o   5   i    l-l    in  to        to  r 

number  (64  bits)  implementation  dictates,  for  efficiency,  the 
use  of  floating  point  hardware  instead  of  the  general  registers 
(32  bits).   Each  w.  can  be  represented  as  a  floating  point  num- 
ber with  leading  zeroes  up  to  the  right  most  hexadecimal  digit 
and  with  an  exponent  of  C,  where  C  is  the  total  number  of  frac- 
tional digits.   (For  the  360/370  long  floating  point  number, 
C  is  14.)   This  scheme,  in  effect,  allows  fixed  point  (integer) 
arithmetric  to  be  performed  with  floating  point  hardware. 

We  now  consider  the  algorithm  to  convert  any  decimal  num- 
ber to  binary  while  retaining  significant  digits.   Let 

E 
N  =  ±  x.,  Xp.  .  .x..  y.,yp.  .  .y,  *  10  .   Clearly,  (j+k)  must  be  less 

or  equal  to  D,  the  maximum  number  of  decimal  digits  permitted. 
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(D  is  15  for  370/360).   The  signs  of  N  and  E  are  saved.   Next, 
the  integers  |E|,  k,  and  x, x? . . .x ,y, Vp . . . yR  are  converted  to 
binary  using  the  previous  algorithm.   If  it  is  possible  to  intro- 
duce a  guard  digit,  then  I,  the  binary  representation  of 

x-.  Xp  .  .  .x  .  y,  y?  .  .  .yR  is  then  adjusted. 

F— k 
The  next  step  is  to  multiply  I  by  the  quantity  10  ~   using 

significant  digit  (SD)  arithmetic.   In  this  way,  the  final  number 

will  have  the  same  significance  as  the  original  number.   This 

step  can  be  separated  into  three  cases,  depending  upon  the 

sign  and  magnitude  of  E.   The  cases  are: 

(1)  E  >  0.   Let  ES  =  I-k.   For  ES  <  0,  divide  I  by  10 ' ES ' - 

For  0  £  ES  <  M,  where  M  is  the  maximum  representable  power  of  ten, 

I  FS  I 

multiply  I  by  10'   1.  (Because  of  unnormalization,  M  is  49  for  an 

imprecise  number  while  M  is  63  for  a  precise  number.) 

I  F  I 

(2)  -M  £  E  <  0.   Divide  I  by  10 '  '  and  then  divide  the 

k 
result  by  10  . 

/  v  M     I E I -M 

(3)  E  <  -M.   Divide  I  successively  by  10  ,  10 ■  '   ,  and 

10  .   In  this  case,  E  has  the  lower  bound  of  smallest  represent- 
able exponent.   (For  the  present  work,  E  must  be  equal  or  greater 
than  -63. ) 

All  the  operations  in  the  above  three  cases  are  performed 
in  significant  digit  arithmetic.   The  last  step  is  to  attach 
the  proper  algebraic  sign  to  the  number. 
4.3   Binary  to  Decimal  Conversion 

Let  g  be  the  number  of  guard  digits  in  the  fractional  part 
of  the  floating  point  number.   (In  the  present  work,  imprecise 
numbers  have  one  guard  digit  and  precise  numbers  have  none.) 
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The  first  step  is  to  convert  the  original  number  to  a  signifi- 
cant integer  multiplied  by  a  power  of  ten.   There  are  three 
cases  depending  on  x,  the  exponent  of  the  number  to  be  conver- 
ted . 

(1)  x  >  (C-g).   Let  ES  =  x  -  (C-g)  and  E(n)  denote  the 
exponent  of  the  normalized  binary  representation  of  10  .   The 
original  number  is  divided  by  10  ,  where  n  is  chosen  such  that 
E(n)  -  1  <  ES  <  E(n+1).   If  the  exponent  of  the  quotient  is 
less  than  (C-g),  then  the  number  is  divided  instead  by  10 
After  conversion  of  the  significant  integer  to  decimal,  the 
result  is  the  decimal  integer  times  10   or  10    . 

(2)  (C-g)  -  G  <  x  <  (C-g),  where  G  is  the  maximum  repre- 
sentable  exponent  in  the  floating  point  number.   The  original 
number  is  multiplied  by  10  ,  where  n  is  chosen  such  that 
E(n-l)  <  ES  <  E(n) .   If  the  exponent  of  the  product  is  less  than 
(C-g) ,  the  number  is  mutliplied  instead  by  10    .   The  result 

is  the  converted  integer  times  10~   or  n0~  ~  . 

M 

(3)  x  <  (C-g)  -  G.   The  number  is  multiplied  by  10  ,  where 

M  is  the  maximum  representable  power  of  ten.   Then  the  resulting 
product  is  treated  as  in  case  (2)  above.   Here  the  result  is  the 
converted  integer  times  10  ^   ,  p  is  the  power  of  ten  used  in 
the  second  multiplication. 

Finally,  if  a  decimal  point  is  to  be  inserted  in  any  de- 
sired position,  the  power  of  ten  is  modified  accordingly. 
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4.4   The  I/O  Routines 

The  input  and  output  of  imprecise  numbers  conform  to  the 
conventions  that  are  used  for  PL/1  E  and  F  formats.   The  special 
symbol  P  is  used  to  distinguish  precise  numbers  from  imprecise 
numbers,  e.g.,  LOOP  and  -1.23P+01.   Like  PL/1,  the  maximum 
number  of  digits  per  item  allowed  is  15,  due  to  the  computer 
word  length.   The  I/O  routines  consist  of  three  modules: 
0L2PILK,  which  accepts  the  input  numbers  as  character  strings 
in  the  PL/1  GET  LIST  manner,  0L2PICP,  which  converts  character 
strings  into  long  floating  number,  and  0L2PIPD,  which  converts 
the  internal  numbers  into  character  strings  for  output.   Both 
0L2PICF  and  0L2PIFD  perform  an  operation  with  a  ten  of  power, 
as  described  by  Kanner ,  in  two  successive  steps  in  the  ten's 
place  and  then  in  the  one's  place  to  preserve  precise  numbers. 

It  was  found  that  a  precise  number  multiplied  by  10   became 

20 

imprecise,  but  the  same  precise  number  multiplied  by  10   and 

10  remained  precise.   This  is  a  result  of  the  computer  word 
length  used. 

Besides  returning  the  correct  significance  for  a  number 
on  output,  0L2PIFD  also  supplies  other  information  to  interpret 
the  number.   There  is  significance  pointer  to  indicate  the  last 
significant  digit  of  the  output  number.   This  is  useful,  for 
example,  when  a  user  does  not  want  to  output  all  significant 
digits  of  a  number.   A  common  user  error  is  to  specify  an  in- 
sufficient field  width  for  output.   In  this  case  an  error  flag 
is  returned  along  with  as  much  of  the  number  as  possible.   For 
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output  of  arrays  it  is  also  useful  to  have  indicators  for  a 
maximum  significant  digit  count  and  a  maximum  field  width. 
This  allows  the  user  to  analyze  his  array  output  with  respect 
to  maximum  significance  and  field  width.   Each  of  these  consid- 
erations are  valuable  and  contribute  to  a  reliable  and  useful 
significant  digit  I/O  package. 
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CHAPTER  5 
DISCUSSION 

Precise  and  imprecise  arithmetic,  as  described  up  to  this 
point,  can  be  used  to  carry  out  more  complicated  mathematical 
operations.   However,  separate  routines  to  handle  operations 
with  complex  numbers  and  matrices  may  be  advisable.   Further  in- 
vestigation is  needed  to  study  the  effects  of  correlated  errors 
with  imprecise  numbers.   More  significance  rules  may  be  required 
or  perhaps,  the  present  significant  rules  must  be  modified. 

5.1  Matrix 

Most  matrix  operations  consist  of  a  series  of  inner  prod- 
ucts of  vectors.   The  inner  product  subroutine  0L2PIVP  was 
written  to  extend  ABC  arithmetic  to  matrix  operations.   Follow- 
ing the  same  approach  of  the  significant  functions,  the  addition 
and  multiplication  are  done  normalized.   At  the  end,  the  result 
is  unnormalized  to  the  correct  number  of  significant  digits  by 
keeping  a  count  of  the  maximum  number  of  leading  zeroes  encount- 
ered.  In  this  way,  the  resulting  inner  product  is  more  accurate 
because  of  less  round-off  errors  with  normalized  numbers. 

5.2  Complex  Numbers 

A  complex  number  can  be  represented  as  a  number  with  real 
and  imaginary  parts.   Routines  using  the  approach  of  the  inner 
product  can  extent  ABC  arithmetic  to  include  complex  numbers. 
Some  definitions  are  necessary:   a  complex  number  is  precise  if 
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and  only  if  both  real  and  imaginary  parts  are  precise;  the  sig- 
nificance of  a  complex  number  is  that  of  the  less  significant 
part;  a  complex  true  zero  consists  of  two  real  true  zeroes. 
The  possibility  of  a  relative  zero  part  and  non-zero  part  must 
be  studied  and  appropriate  definitions  implemented. 
5.3   Some  Problems  with  Significant  Digit  Arithmetic 

The  adjustment  rules  for  imprecise  numbers  are  based  on 
the  assumption  that  their  associated  errors  are  statistically 
independent.   If  their  errors  are  correlated,  the  resulting 
imprecise  number  may  have  a  representation  which  is  not  a  re- 
liable measure  of  significance.   An  example  of  correlated  error 
is  found  in  the  multiplication  of  an  imprecise  number  by  itself. 
The  resulting  product  should  be  more  significant  since  the 
approximation  error  was  squared.   The  auxiliary  rule  for  multi- 
plication attempts  to  correct  this  situation  by  increasing  the 
number  of  significant  digits  by  one.   Similarly,  repeated  sum- 
mation of  an  imprecise  number  by  itself,  under  the  present  ad- 
justment rule,  would  allow  the  number  of  significant  digits  to 
increase.   A  combination  of  these  correlated  errors  can  be  found 
in  the  vector  norm  operation.   Thus,  correlated  errors  in  signi- 
ficant digit  arithmetic  should  be  studied. 

To  complete  the  implementation  of  significant  digit  arith- 
metic, routines  to  compare  two  operands  are  needed.   Detection 
of  a  relative  and  true  zero  is  done  by  check  for  zero  fraction. 
If  the  exponent  is  also  all  zero,  then  the  number  is  a  true 
zero.   To  compare  two  relative  zeroes,  the  exponents  are  examined. 
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Since  relative  zeroes  have  no  significant  digits,  only  the  ex- 
ponents of  a  non-zero  and  a  relative  zero  can  be  compared.   True 
zero  would  always  be  "lower"  than  any  number  which  is  not  a  true 
zero.   Care  must  be  taken  to  avoid  using  the  guard  bits  of  an 
imprecise  number  in  a  compare.   After  removing  the  guard  bits 
and  normalizing,  comparison  of  imprecise  numbers  can  be  made. 
Finally,  precise  numbers  are  compared  with  any  special  treatment. 

As  noted  in  Chapter  3,  the  present  implementation  of  func- 
tional values  for  relative  zero  arguments  did  not  account  for 
the  possibility  of  repeated  function  evaluations  using  relative 
zeroes.   Each  function  should  be  studied  separately  to  provide 
a  more  comprehensive  treatment  of  relative  zero  arguments  and 
correlated  errors  in  functional  evaluations  in  significant  digit 
arithmetic . 

There  are  other  methods  of  error  monitoring.   The  interval 
arithmetic  of  Moore  (19)  also  provides  error  bounds  on  computed 
results.   Metropolis  (17)  gives  techniques  of  analyzing  algorithms 
in  unnormalized  arithmetic.   Computational  errors  are  also  dis- 
cussed in  (21).   Metropolis   (  9  )  states  "...  different  algor- 
ithms for  the  same  problem  give  different  representations  is 
connected  with  the  observation  that  mathematically  equivalent 
approaches  are  not  computationally  equivalent."   Significant 
digit  arithmetic  is  one  way  to  aid  the  error  analysis  of  a  com- 
putational algorithm.   It  is  hoped  that  present  implementation 
will  aid  further  study  and  experimentation  with  imprecise  numbers 
and  significant  digit  arithmetic. 
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APPENDIX  A 
EXAMPLES  OF  SIGNIFICANT  DIGIT  ARITHMETIC 
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APPENDIX  A 

This  appendix  contains  several  examples  of  significant 
digit  arithmetic .   Each  program  illustrates  how  a  software  sys- 
tem may  generate  the  necessary  linkage  or  CALLS  to  the  basic 
significant  digit  routines.   The  unnormalized  significance  out- 
put was  provided  by  0L2PIPD.   The  significance  pointer  is  the 
position  of  the  last  significant  digit,  counting  from  the  left- 
most position  of  the  field.   This  significance  pointer  can  also 
return  negative  numbers:   -3  for  true  and  relative  zeroes,  -2 
for  insufficient  field  width  (severe  error),  and  -1  to  indicate 
that  not  all  of  the  significant  digits  in  the  fraction  were  dis- 
played for  the  given  field  width  (warning  error) .   The  informa- 
tion in  the  significance  pointer  may  be  incorporated  in  a  more 
elaborate  I/O  routine  which  can  place  a  marker  such  as  the 
underscore  character  under  the  last  significant  digit  of  an 
output  item. 
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APPENDIX  B 
EXAMPLES  OF  DOCUMENTATION  OF  THE  BASIC  ROUTINES 
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APPENDIX  B 

This  appendix  contains  the  assembly  listings  of  several 
significant  digit  routines:  @0L2ADD,  @0L2MPY,  and  0L2PIVP. 
One  of  the  assembly  language  macros  used  to  generate  the  sig- 
nigicant  digit  function  is  also  included.  It  should  be  noted 
that  in  its  present  form,  0L2PIVP  (inner  product)  assumes  the 
elements  in  the  input  vector  occupying  consecutive  locations. 
The  indexing  scheme  can  be  modified  to  handle  more  complicated 
storage  schemes  as  found  in  the  array  structures  of  OL/2  lan- 
guage (20)  . 
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